;;------------------------------------------------------------------------- 
;;   - Drawing box plots
;;   - Using a box plot to show the the median, minimum/maximum value, and the 25th/75th percentiles of two timeseries
;;   - Adding text to a box plot
;;   - Sorting data
;;   - Setting the mininum/maximum value of the Y axis in a box plot
;;   - Generating dummy data using "random_normal"
;;------------------------------------------------------------------------- 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$KAILIBROOT/lib_radon.ncl"
load "$KAILIBROOT/lib_boxplot.ncl"

begin

  cpp    = "" 

  ss = (/ "ALL","DJF","JJA","MAM","SON" /)

  sname = asciiread("radon_sitename.dat",-1,"string")

  ns = dimsizes(sname)

  plot = new(ns,graphic)  
 
;;------------------------------------------------------------------------- 
;; load settings by user 
;;------------------------------------------------------------------------- 

;;------------------------------------------------------------------------- 
;; filename 
;;------------------------------------------------------------------------- 
  filenama = "box_data_ANN.nc" 
  filenamb = "box_data_DJF.nc" 
  filenamc = "box_data_MAM.nc"
  filenamd = "box_data_JJA.nc"
  filename = "box_data_SON.nc"

  haveobs = isfilepresent(filenama) 

  fla = addfile(filenama,"r") 
  flb = addfile(filenamb,"r") 
  flc = addfile(filenamc,"r") 
  fld = addfile(filenamd,"r") 
  fle = addfile(filename,"r") 

;;------------------------------------------------------------------------- 
;; read data 
;;------------------------------------------------------------------------- 
  ppa = fla->boxa
  ppb = flb->boxa
  ppc = flc->boxa
  ppd = fld->boxa
  ppe = fle->boxa

  paa = fla->boxb
  pab = flb->boxb
  pac = flc->boxb
  pad = fld->boxb
  pae = fle->boxb

  is = 0 

  print("") 
  print("read data") 
  print("") 


  do is = 0,ns-1 

;;------------------------------------------------------------------------- 
;; plot setting 
;;------------------------------------------------------------------------- 
  wks = gsn_open_wks("pdf","plotbox_"+sname(is)+cpp)                ; create postscript file

  res            = True                         ; plot mods desired
 ;res@tmXBLabels = (/"OBS","E0","E1","E2","E3"/) ; labels for each box
  res@tmXBLabels = (/"OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3", \ 
                     "OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3", \
                     "OBS","E0","E2","E3"/) ; labels for each box
  res@tmXBLabels = (/"","    ANN","","", \
                     "","    DJF","","", \
                     "","    MAM","","", \
                     "","    JJA","","", \
                     "","    SON","",""  /)

  ;;res@tiMainString = sname(is) ;;+ ", "+season
  ;;res@gsnLeftString = sname(is) ;;+ ", "+season

 ;res@tiXAxisString    = ""
  res@tiYAxisString    = "Radon activity (mBq m~S~-3~N~STP) "    ; y-axis label     
  res@tiYAxisString    = "IPR (cm~S~-3~N~s~S~-1~N~) "    ; y-axis label     

  ;;res@trYMinF = -900.
  res@trYMinF = 0.



;;------------------------------------------------------------------------- 
;; merge data  
;;------------------------------------------------------------------------- 

  nep = 4 
  npt = 5 
  nsn = 5 

  ;ncl0 = 121 ;
  ;ncl1 = 4 ;
  ;ncl2 = 5 ;

  xpp = new((/ns,nep,nsn,npt/),"float") 
  xpa = new((/ns,nep,nsn/),"float") 

  nsn = 5 
 
  do isn = 0, nsn-1
     xpp(:,:,0,:) = ppa(:,:,:)  
     xpp(:,:,1,:) = ppb(:,:,:)  
     xpp(:,:,2,:) = ppc(:,:,:)  
     xpp(:,:,3,:) = ppd(:,:,:)  
     xpp(:,:,4,:) = ppe(:,:,:)  
  end do 
 
  xpa(:,:,0) = paa(:,:)  
  xpa(:,:,1) = pab(:,:)  
  xpa(:,:,2) = pac(:,:)  
  xpa(:,:,3) = pad(:,:)  
  xpa(:,:,4) = pae(:,:)  

  pp = new((/ns,nep*nsn,npt/),"float") 
  pa = new((/ns,nep*nsn/),"float") 

printVarSummary(pp) 
printVarSummary(xpp) 

  isp = 0 

  do isn = 0, nsn-1 ;; season second 
  do iep = 0, nep-1 ;; exp first  
     print("isn,iep,isp: "+isn+" "+iep+" "+isp) 
     ;;     isn,iep,isp: 4 3 20
     pp(:,isp,:) = xpp(:,iep,isn,:) 
     pa(:,isp)   = xpa(:,iep,isn) 
     isp = isp + 1  
  end do 
  end do 

  ;pp = RN222_ionprod(pp)  
  ;pa = RN222_ionprod(pa)  

  pp = where(pp.lt.1.E-06,-10000,pp) 
  pa = where(pa.lt.1.E-06,-10000,pa) 


 ;ppx = ispan(0,nep*nsn-1,1) 
  ppx = new( nep*nsn, float )
  do isn=0,nsn-1
  do iep=0,nep-1
      ppx(isn*nep+iep) = 5 + isn*(nsn*0.6) + iep*0.4
  end do
  end do

  print("") 
  print("compute vars") 
  print("") 

  ismax = max(pp(is,:,:)) .gt. max(pa(is,:)) 

  res@trXMinF = 2  
  res@trXMaxF = 22  

  if(ismax) then  
  res@trYMaxF = max(pp(is,:,:)) * 1.2 
  else
  res@trYMaxF = max(pa(is,:)) * 1.2 
  end if 

  res@mpShapeMode  = "FreeAspect"
  res@vpWidthF      = 0.8
  res@vpHeightF     = 0.4


  ;if(sname(is).eq."grifton496".and.season.eq."JF") then 
  ;   res@trYMaxF =  8000.
  ;end if 
  ;if(sname(is).eq."grifton496".and.season.eq."MAM") then 
  ;   res@trYMaxF =  5000.
  ;end if 

  if (False) then 

  if(sname(is).eq."grifton") then 
     res@trYMaxF =  10000.
  end if 
  if(sname(is).eq."capepoint") then 
     res@trYMaxF =  4000.
  end if 
  if(sname(is).eq."capegrim") then 
     res@trYMaxF =  4000.
  end if 
  if(sname(is).eq."hongkong") then 
     res@trYMaxF =  10000.
  end if 
  if(sname(is).eq."gosan") then 
     res@trYMaxF =  5000.
  end if 
  if(sname(is).eq."freiburg") then 
     res@trYMaxF =  15000.
  end if 
  if(sname(is).eq."schauinsland") then 
     res@trYMaxF =  5000.
  end if 
  if(sname(is).eq."bermuda") then 
     res@trYMaxF =  2000.
     ;res@trYMinF = -400.
  end if 

  ;if(sname(is).eq."bermuda".and.season.eq."DJF") then 
  ;   res@trYMaxF =  2500.
  ;   res@trYMinF = -400.
  ;end if 
  ;if(sname(is).eq."bermuda".and.season.eq."JJA") then 
  ;   res@trYMaxF =  1500.
  ;   res@trYMinF = -200.
  ;end if 


  end if 


;**********************************************
; legend 
;**********************************************
 ;res@pmLegendDisplayMode    = "Always"              ; turn on legend
 ;res@pmLegendSide           = "Top"                 ; Change location of 
 ;res@pmLegendParallelPosF   = .90                   ; move units right
 ;res@pmLegendOrthogonalPosF = -0.8                  ; more neg = down
 ;res@pmLegendWidthF         = 0.12                  ; Change width and
 ;res@pmLegendHeightF        = 0.25                  ; height of legend.
 ;res@lgLabelFontHeightF     = .02                   ; change font height
 ;res@lgPerimOn              = False                 ; no box around
 ;res@xyExplicitLegendLabels = (/"20N","30N","40N","50N"/)


;**********************************************
; resources for polylines that draws the boxes
;**********************************************  
  llres                   = True                        
  llres@gsLineThicknessF  = 2.5                 ; line thickness 
;**********************************************
; resources that control color and width of boxes
;**********************************************  
  opti          = True                  
  opti@boxWidth = .25                           ; Width of box (x units)
  opti@boxColors = (/"black","blue","red","green", \  
                     "black","blue","red","green", \
                     "black","blue","red","green", \ 
                     "black","blue","red","green", \
                     "black","blue","red","green"  /)     ; Color of box(es)

  print("") 
  print("making plot") 
  print("") 

;**********************************************
; make plot 
;**********************************************
  bv = pp(is,:,:) 

  printVarSummary(ppx) 
  printVarSummary(bv) 

  do i = 0,nsn*nep-1  
     print(ppx(i)+"   "+bv(i,:)) 
  end do 

  res@tmXBMajorOutwardLengthF = 0. ;;0.021   ; set tickmarks 
  res@tmXBMinorOutwardLengthF = 0. ;;0.015
  res@tmXBMajorLengthF = 0.006
  res@tmXBMinorLengthF = 0.006
  res@tmYLMajorLengthF = 0.010
  res@tmYLMinorLengthF = 0.006

  plot(is) = boxplot2(wks,ppx,bv,opti,res,llres)      

  print("") 
  print("plot defined") 
  print("") 


;**********************************************
; adding additional text to the plot
;**********************************************  
;  tres = True
;  str  = "E0: WCRP1995 ~C~"+\
;         "E2: Schery (1998) scaled by 1.6 ~C~"+ \
;         "E3: Merged radon flux map" 
;         ;Szegvary et al. (2007)
;         ;Zhuo et al. (2008)
;         ;Griffiths et al. (2010)
;         ;Schery (1998)
;  tres@txFontHeightF = 0.015
;  gsn_text_ndc(wks,str,0.60,0.73,tres)

;;------------------------------------------------------------------------- 
;; add some polymarkers
;;------------------------------------------------------------------------- 
  mres               = True                     ; marker mods desired
  mres@gsMarkerIndex = 1                        ; polymarker style
  mres@gsMarkerSizeF = 40.                      ; polymarker size
  mres@gsMarkerColor = "black"                  ; polymarker color

  print("") 
  print("making markers") 
  print("") 
  dum = new((/ns,nsn*nep/),graphic) 

  ;ppx = new( nep*nsn, float )
  ;do isn=0,nsn-1
  ;do iep=0,nep-1
  ;    ppx(isn*nep+iep) = 5 + isn*(nsn*0.6) + iep*0.4
  ;end do
  ;end do

  do isp = 0,nsn*nep-1 
     dum(is,isp) = gsn_add_polymarker(wks,plot(is),ppx(isp),pa(is,isp),mres) 
  end do 

  draw(wks)                                    ; boxplot does not call these
  frame(wks)                                   ; for you

  delete(wks) 

  end do ;; ns 


end






